Load data

Intensity and shapefile data is handled here.

# make dates
info = gdalinfo("VV_2017_clip.tif")
descr = info[grepl(info, pattern = "Description = ")]
descr = gsub(descr, pattern = "Description = ", replacement = "")
dates2017 <- as_date(descr)

info = gdalinfo("VV_clip.tif")
descr = info[grepl(info, pattern = "Description = ")]
descr = gsub(descr, pattern = "Description = ", replacement = "")
dates2018 <- as_date(descr)
# function takes a shape as an AOI and a year (2017 or 2018)
loadPolygon <- function(shape, year) {
  if(year == 2017) {
    int_VV <- int_VV_2017
    int_VH <- int_VH_2017
    dates <- dates2017
  } else {
    if(year == 2018) {
      int_VV <- int_VV_2018
      int_VH <- int_VH_2018
      dates <- dates2018
    } else {
      return("wrong year entered!")
    }
  }
  int_vv <- st_as_stars(int_VV[shape])
  int_vh <- st_as_stars(int_VH[shape])
  int_vv <- st_set_dimensions(int_vv, 3, val = dates, names = "time")
  int_vh <- st_set_dimensions(int_vh, 3, val = dates, names = "time")
  # c()
  comb <- c(int_vv, int_vh, along="bands")
  # switch bands to attributes, this is more of a design choice.
  comb_split <- split(comb, "bands")
  names(comb_split) <- c("VV", "VH")
  return(comb_split)
}

Combined Water Threshold Training

Detection Goal

area1 <- loadPolygon(shape[124,], 2018)
area2 <- loadPolygon(shape[124,], 2017)
# this kind of water puddle:
plot(area1[1,,,1])

Make Data From Both Years

Before, thresholds were calculates for single years, resulting in -15.3 for 2018 and about -18 for 2017.

# prepare data
val.svm$type <- as.factor(val.svm$type)
attach(val.svm)
x <- subset(val.svm, select=c(-type, -sitch))
y <- type
# make model
svmmod <- svm(x,y)
# summary(svmmod)
# make prediction, confusion matrix
pred <- predict(svmmod,x)
table(pred, y)
##        y
## pred    field water
##   field  4019    68
##   water     1  1948
# threshold is
svmmod$x.scale$`scaled:center`
## [1] -17.62537
# plot classified with threshold calculated by SVM
plot(area1[1,,,1], col = c("white", "black"), breaks = c(-35, svmmod$x.scale$`scaled:center`, 10))
# vs old threshold
# plot(all[1,,,1], col = c("white", "black"), breaks = c(-35, thresh$threshold, 10))
ggplot(val.svm, aes(x=sitch, y=VV)) +
  geom_boxplot() +
  geom_hline(yintercept = svmmod$x.scale$`scaled:center`)

Look at Training Data

## The following objects are masked from val.svm (pos = 3):
## 
##     sitch, type, VV
##        y
## pred    field water
##   field  4060   130
##   water     5  2000
## [1] -17.55717

Create New Threshold - Sum Raster

Plot, Threshold -17

allRas2017 <- raster("./minus17/allRas2017.tif")
allRas2018 <- raster("./minus17/allRas2018.tif")
colo = viridisLite::inferno(29)
breaks = seq(0, 29, 1)
image(allRas2017, col = colo, breaks = breaks, main="2017")
colo = viridisLite::inferno(31)
breaks = seq(0, 31, 1)
image(allRas2018, col = colo, breaks = breaks, main="2018")

Comparison to Threshold = -15.3

allRas2017 <- raster("./allRas2017.tif")
allRas2018 <- raster("./allRas2018.tif")
colo = viridisLite::inferno(29)
breaks = seq(0, 29, 1)
image(allRas2017, col = colo, breaks = breaks, main="2017")
colo = viridisLite::inferno(31)
breaks = seq(0, 31, 1)
image(allRas2018, col = colo, breaks = breaks, main="2018")

Zoom, Threshold -17

allRas2017 <- read_stars("./minus17/allRas2017.tif")
allRas2018 <- read_stars("./minus17/allRas2018.tif")
st_crs(allRas2017) <- "EPSG:25832"
st_crs(allRas2018) <- "EPSG:25832"
colo = viridisLite::inferno(29)
breaks = seq(0, 29, 1)
image(allRas2017[study_area[1,]], col = colo, breaks = breaks, main="2017")
colo = viridisLite::inferno(31)
breaks = seq(0, 31, 1)
image(allRas2018[study_area[1,]], col = colo, breaks = breaks, main="2018")
st_crs(allRas2018) <- "EPSG:25832"

Load RainData

## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

Monitor Moorflächen

# cut shapefile to area of interest
inter_shape <- shape[lengths(st_intersects(shape, study_area[1,])) != 0,]
plot(inter_shape)

# load
moor2017 <- loadPolygon(inter_shape, 2017)
plot(moor2017)

# set - 23 as Threshold for open water
moor2017[moor2017[1,,,] < -23.5] <- NA
# plot(moor2017)

makeTimeSeries <- function(stars, shape) {
  # for each shape
  for (i in 1:length(shape$Förderung)) {
    cut <- stars[shape[i,]]
    cut <- as.data.frame(st_apply(cut, "time", mean, na.rm = TRUE))
    
    if(!is.na(shape[i,]$Förderung)) {
      names(cut) <- c("time", paste0("F", i, "VV"), paste0("F", i, "VH"))
    } else {
      names(cut) <- c("time", paste0(i, "VV"), paste0(i, "VH")) 
    }
    
    # at first iteration
    if(i == 1) {
      # create master df
      df <- cut
    # at every ther iteration
    } else {
      # append
      df <- cbind(df, cut[,2:3])
    }
  }
  return(df)
}

plotTimeSeries <- function(df, year) {
  if(year==2017) {
    rain.select <- rain.all.2017.df
  } else {
    if(year==2018) {
      rain.select <- rain.all.2018.df
    }
    else {
      if(class(year) == "data.frame") {
        rain.select <- year
      }
    }
  }
  g <- ggplot() + geom_bar(aes(x=rain.select$time, y=rain.select$mean), stat='identity', fill = "red", alpha = 0.8) + coord_cartesian(ylim=c(0,35)) + ggtitle("Mean of Moorflächen (*2) plotted against Precipitation (3 day sum) \nblue = in Förderung, black = nicht in Förderung") + xlab("Time") + ylab("Precipitation in mm, Backscatter VV")
  #+ scale_y_continuous(sec.axis = sec_axis(~. /10 + 22, name = "Intensity + 22"))
  # select VV columns
  for (i in seq(2, ncol(df), 2)) {
    if(names(df)[i] == paste0("F", i/2, "VV")) {
      g <- g + geom_line(aes_string(x=df[,1], y = df[,i]*2 + 44, alpha = 0.7), color = "blue")
      g <- g + geom_point(aes_string(x=df[,1], y = df[,i]*2 + 44, alpha = 0.7), color = "blue")
    } else {
      g <- g + geom_line(aes_string(x=df[,1], y = df[,i]*2 + 44), alpha = 0.3)
    }
  }
  print(g)
}

df <- makeTimeSeries(moor2017, inter_shape)
plotTimeSeries(df, 2017)

# load
moor2018 <- loadPolygon(inter_shape, 2018)

# set - 23 as Threshold for open water
moor2018[moor2018[1,,,] < -23.5] <- NA

plotTimeSeries(makeTimeSeries(moor2018, inter_shape), 2018)

dff <- rbind(makeTimeSeries(moor2017, inter_shape), makeTimeSeries(moor2018, inter_shape))
plotTimeSeries(dff, year = rbind(rain.all.2017.df, rain.all.2018.df))

dfff <- dff[,c(1,seq(2,ncol(dff), 2))]
gef <- c(1)
ngef <- c(1)
for (i in 2:length(names(dfff))) {
  str <- paste0("F", i - 1, "VV")
  if(names(dfff)[i] == str) {
    gef <- c(gef, i)
  } else {
    ngef <- c(ngef, i)
  }
}

inF <- dfff[,gef]
ninF <- dfff[,ngef]

inF$Mean <- rowMeans(inF[,2:ncol(inF)])
ninF$Mean <- rowMeans(ninF[,2:ncol(ninF)])

dfM <- cbind(inF[,c(1, ncol(inF))], ninF$Mean, ninF$Mean)
names(dfM) <- c("time", "F1VV", "1XX", "1VV")

plotTimeSeries(dfM, year = rbind(rain.all.2017.df, rain.all.2018.df))
## Warning in if (year == 2017) {: the condition has length > 1 and only the first
## element will be used
## Warning in if (year == 2018) {: the condition has length > 1 and only the first
## element will be used

Load New Rain Data

Load RainData

## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

Look at Means

head(dfM)
##         time      F1VV       1XX       1VV
## 1 2017-01-08 -14.97066 -15.63378 -15.63378
## 2 2017-02-01 -12.78420 -12.97320 -12.97320
## 3 2017-02-13 -15.60427 -14.87243 -14.87243
## 4 2017-02-25 -11.78134 -12.13760 -12.13760
## 5 2017-03-09 -11.77774 -12.27017 -12.27017
## 6 2017-03-21 -11.97166 -11.94932 -11.94932
# delete the column needed for plotting function
dfM <- dfM[,c(1,2,4)]
# all acquisitions are 12 days apart, the first one is 8.1.2017

# rain
rain.all <- rbind(rain.all.2017.df, rain.all.2018.df)

# make sum rain vector
vec <- c(rain.all[1,2] + rain.all[2,2] + rain.all[3,2] + rain.all[4,2] + rain.all[5,2] + rain.all[6,2] + rain.all[7,2] + rain.all[8,2])
for (i in 1:nrow(dfM)-1) {
  offset <- 8
  inter <- 12
  count <- i - 1
  start <- count*inter + offset + 1
  end <- i * inter + offset
  sum <- 0
  # print(c(start, end))
  for (j in start:end) {
    sum <- sum + rain.all[j,2]
  }
  vec <- c(vec, sum)
}

df <- cbind(dfM, vec)
names(df) <- c("time", "FVV", "nFVV", "prec")

ggplot(df, aes(x = time)) + 
  geom_bar(aes(x = time - 6, y = prec), stat = 'identity', fill = "lightblue", alpha = 0.8) + 
  geom_line(aes(y = FVV * 4 + 70, color = "in Förderung")) +
  geom_point(aes(y = FVV * 4 + 70, color = "in Förderung")) +
  geom_line(aes(y = nFVV * 4 + 70, color = "nicht in Förderung")) +
  coord_cartesian(ylim = c(0,60)) + 
  ggtitle("Mean of Moorflächen and Precipitation") + xlab("Time") +
  ylab("Precipitation in mm/m²") + 
  scale_y_continuous(sec.axis = sec_axis(~. *0.25 -17.5, name = "Intensity in dB")) + 
  scale_color_manual(name = "Förderung", values = c("in Förderung"="blue", "nicht in Förderung"="black")) +
  theme(legend.position = "bottom")

# make sum rain vector
vec <- c(rain.all[1,2] + rain.all[2,2] + rain.all[3,2] + rain.all[4,2] + rain.all[5,2] + rain.all[6,2] + rain.all[7,2] + rain.all[8,2])
for (i in 1:nrow(dfM)-1) {
  offset <- 8
  inter <- 12
  count <- i - 1
  start <- count*inter + offset + 1
  end <- i * inter + offset
  sum <- 0
  # print(c(start, end))
  weight <- 1
  for (j in start:end) {
    sum <- sum + (1/weight * rain.all[j,2])
    weight <- weight + 1
  }
  vec <- c(vec, sum)
}

gf <- cbind(dfM, vec)
names(gf) <- c("time", "FVV", "nFVV", "prec")

ggplot(gf, aes(x = time)) + 
  geom_bar(aes(x = time - 6, y = prec), stat = 'identity', fill = "lightblue", alpha = 0.8) + 
  geom_line(aes(y = FVV * 1.5 + 28, color = "in Förderung")) +
  geom_point(aes(y = FVV * 1.5 + 28, color = "in Förderung")) +
  geom_line(aes(y = nFVV * 1.5 + 28, color = "nicht in Förderung")) +
  coord_cartesian(ylim = c(0,25)) + 
  ggtitle("Mean of Moorflächen and Precipitation, Weight: 1 / # of Day") + xlab("Time") +
  ylab("Precipitation in mm/m²") + 
  scale_y_continuous(sec.axis = sec_axis(~. *2/3 - 56/3, name = "Intensity in dB")) + 
  scale_color_manual(name = "Förderung", values = c("in Förderung"="blue", "nicht in Förderung"="black")) +
  theme(legend.position = "bottom")

Normalized Difference VV - Rain

head(df)
##         time       FVV      nFVV       prec
## 1 2017-01-08 -14.97066 -15.63378 15.3018018
## 2 2017-02-01 -12.78420 -12.97320 10.7297298
## 3 2017-02-13 -15.60427 -14.87243  5.3207207
## 4 2017-02-25 -11.78134 -12.13760  0.8630631
## 5 2017-03-09 -11.77774 -12.27017 22.7909910
## 6 2017-03-21 -11.97166 -11.94932 22.4567567
dff <- df
dff$FVVbyprec <- (25 + dff$FVV - dff$prec) / (25 + dff$FVV + dff$prec)


ggplot(dff, aes(x = time)) + 
  geom_line(aes(x = time - 6, y = FVVbyprec)) +
  
  geom_bar(aes(x =time - 6, y = prec), stat = 'identity', fill = "lightblue", alpha = 0.8) + 
  geom_line(aes(y = FVV * 4 + 70, color = "in Förderung")) +
  geom_point(aes(y = FVV * 4 + 70, color = "in Förderung")) +
  geom_line(aes(y = nFVV * 4 + 70, color = "nicht in Förderung")) +
  coord_cartesian(ylim = c(-5,40)) + 
  ggtitle("Mean of Moorflächen and Precipitation") + xlab("Time") +
  ylab("Precipitation in mm/m²") + 
  scale_y_continuous(sec.axis = sec_axis(~. *0.25 -17.5, name = "Intensity in dB")) + 
  scale_color_manual(name = "Förderung", values = c("in Förderung"="blue", "nicht in Förderung"="black")) +
  theme(legend.position = "bottom")

dff$Fchange <- 1:nrow(dff)
for (i in 2:nrow(dff)) {
  dff[i,6] <- (dff[i,2] - dff[i-1,2]) / (dff[i,2] + dff[i-1,2])
}
dff$Rchange <- 1:nrow(dff)
for (i in 2:nrow(dff)) {
  dff[i,7] <- (dff[i,4] - dff[i-1,4]) / (dff[i,4] + dff[i-1,4])
}
dff$FRchange <- 1:nrow(dff)
for (i in 2:nrow(dff)) {
  dff[i,8] <- (dff[i,6] - dff[i-1,7]) / (dff[i,6] + dff[i-1,7])
}

# normalized change of normalized change over time with previous time step

ggplot(dff, aes(x = time)) + 
  geom_line(aes(x = time - 6, y = FRchange)) +
  
  geom_bar(aes(x =time - 6, y = prec), stat = 'identity', fill = "lightblue", alpha = 0.8) + 
  geom_line(aes(y = FVV * 4 + 70, color = "in Förderung")) +
  geom_point(aes(y = FVV * 4 + 70, color = "in Förderung")) +
  geom_line(aes(y = nFVV * 4 + 70, color = "nicht in Förderung")) +
  coord_cartesian(ylim = c(-5,40)) + 
  ggtitle("Mean of Moorflächen and Precipitation") + xlab("Time") +
  ylab("Precipitation in mm/m²") + 
  scale_y_continuous(sec.axis = sec_axis(~. *0.25 -17.5, name = "Intensity in dB")) + 
  scale_color_manual(name = "Förderung", values = c("in Förderung"="blue", "nicht in Förderung"="black")) +
  theme(legend.position = "bottom")

cumulative sum

# make cumsum
dff$Rsum <- cumsum(dff$prec - mean(dff$prec[1:30]))

# start from mean prec
dff$prec_1 <- c(mean(dff$prec[1:30]), dff$prec[2:60])
dff$Rsum_1 <- cumsum(dff$prec_1 - mean(dff$prec[1:30]))
# calc means
mean(dff$prec) # all
## [1] 19.53345
mean(dff$prec[1:30]) # 2017
## [1] 26.89676
mean(dff$prec[31:60]) # 2018
## [1] 12.17015
dff$Rsum <- dff$Rsum - mean(dff$prec)

# ggplot(dff[1:30,], aes(x = time)) +
#   geom_bar(aes(y = Rsum / 50), stat = 'identity') +
#   geom_line(aes(y = FVV + 11.5)) +
#   ggtitle("cumsum rain, scaled")

# ggplot(dff[1:30,], aes(x = time)) +
#   geom_line(aes(y = Rsum / 50), color="lightblue") +
#   geom_line(aes(y = FVV + 11.5)) +
#   ggtitle("cumsum rain, scaled")

# R sum 1
ggplot(dff[1:30,], aes(x = time)) +
  geom_bar(aes(y = Rsum_1 / 50), stat = 'identity') +
  geom_line(aes(y = FVV + 12)) +
  ggtitle("cumsum rain, scaled, starting at 1 = mean")

ggplot(dff[1:30,], aes(x = time)) +
  geom_line(aes(y = Rsum_1 / 50), color="lightblue") +
  geom_line(aes(y = FVV + 12)) +
  ggtitle("cumsum rain, scaled, starting at 1 = mean")

# weighted stuff
gf$Rsum <- cumsum(gf$prec - mean(gf$prec[1:30]))

ggplot(gf[1:30,], aes(x = time)) +
  geom_bar(aes(y = Rsum / 10), stat = 'identity') +
  geom_line(aes(y = FVV + 11.5)) +
  ggtitle("cumsum weighted rain, 2017")

ggplot(gf[1:30,], aes(x = time)) +
  geom_line(aes(y = Rsum / 10), color="lightblue") +
  geom_line(aes(y = FVV + 11.5)) +
  ggtitle("cumsum weighted rain, 2017")

# 2018
dff$Rsum[31:60] <- cumsum(dff$prec[31:60] - mean(dff$prec[31:60]))

ggplot(dff[31:60,], aes(x = time)) +
  geom_line(aes(y = Rsum / 4), color="lightblue") +
  geom_line(aes(y = FVV + 17)) +
  ggtitle("cumsum rain 2018, scaled, - mean 2018")

ggplot(dff[1:30,], aes(x = time)) + 
  geom_bar(aes(x =time - 6, y = prec), stat = 'identity', fill = "lightblue", alpha = 0.8) + 
  geom_line(aes(y = FVV * 4 + 70, color = "in Förderung")) +
  geom_point(aes(y = FVV * 4 + 70, color = "in Förderung")) +
  geom_line(aes(y = nFVV * 4 + 70, color = "nicht in Förderung")) +
  geom_hline(yintercept = mean(dff$prec[1:30]), color = "lightblue") +
  coord_cartesian(ylim = c(-5,40)) + 
  ggtitle("Mean of Moorflächen and Precipitation") + xlab("Time") +
  ylab("Precipitation in mm/m²") + 
  scale_y_continuous(sec.axis = sec_axis(~. *0.25 -17.5, name = "Intensity in dB")) + 
  scale_color_manual(name = "Legend", values = c("in Förderung"="blue", "nicht in Förderung"="black", "Mean of Precipitation 2017"="lightblue")) +
  theme(legend.position = "bottom")

# dff$prec_2 <- c(mean(dff$prec), dff$prec[2:60])
dff$Rsum_2 <- cumsum(dff$prec - 23)

ggplot(dff, aes(x = time)) +
  geom_bar(aes(y = Rsum_2 / 50), stat='identity') +
  geom_line(aes(y = FVV * 1.2 + 15))

statistics

acf(dff$prec)
pacf(dff$prec)
acf(dff$FVV)
pacf(dff$FVV)

cor(dff$FVV, dff$prec)
## [1] 0.256756
ccf(dff$FVV, dff$prec)

R_sum1 is the cumulated sum minus the prec mean of 2017, but with pretending that the first precipitation was equal to mean (as a “start” value)

ccf(dff$FVV, dff$Rsum)
ccf(dff$FVV, dff$Rsum_1)
acf(dff$Rsum_1)

dff$acq <- 1:nrow(dff)
f <- function(x) sum((dff$FVV - (x[1] + x[2] * sin(pi * (dff$acq+x[3])/30)))^2)
nlm(f,c(0,0,0))
## $minimum
## [1] 82.44628
## 
## $estimate
## [1] -13.386189   1.610426  -8.671524
## 
## $gradient
## [1] -1.111501e-06 -1.299817e-05  7.397523e-06
## 
## $code
## [1] 1
## 
## $iterations
## [1] 12
dff$FVV.per <- -13.386 + 1.61 * sin(pi*(dff$acq - 8.6715)/30)

ggplot(dff, aes(x = time)) +
  geom_line(aes(y = FVV)) + 
  geom_line(aes(y = FVV.per, color = "red"))

ggplot(dff, aes(x = time)) +
  geom_bar(aes(y = prec / 10), stat = 'identity') + 
  geom_line(aes(y = FVV - FVV.per + 1.8))

an <- dff$FVV - dff$FVV.per

an.ar5 = arima(an, c(5,0,0))
an.ar5
## 
## Call:
## arima(x = an, order = c(5, 0, 0))
## 
## Coefficients:
##          ar1     ar2      ar3      ar4      ar5  intercept
##       0.3567  0.2003  -0.2136  -0.0615  -0.0723     0.0182
## s.e.  0.1295  0.1439   0.1624   0.1626   0.1580     0.1705
## 
## sigma^2 estimated as 1.056:  log likelihood = -87,  aic = 187.99
acf(an, type="partial")

arima(an, c(1,0,0))$aic
## [1] 185.1914
dfh <- dff[1:30,]
n <- 4
f <- function(x) sum((dfh$FVV - (x[1] + x[2] * sin(pi * (dfh$acq+x[3])/n)))^2)
mod <- nlm(f,c(0,0,0))
dfh$FVV.per <- mod$estimate[1] + mod$estimate[2] * sin(pi*(dfh$acq + mod$estimate[3])/n)

 ggplot(dfh, aes(x = time)) +
   geom_line(aes(y = FVV)) + 
   geom_line(aes(y = FVV.per, color = "red"))

n <- 30
f <- function(x) sum((dff$prec - (x[1] + x[2] * sin(pi * (dff$acq+x[3])/n)))^2)
mod <- nlm(f,c(0,0,0))
dff$prec.per <- mod$estimate[1] + mod$estimate[2] * sin(pi*(dff$acq + mod$estimate[3])/n)

ggplot(dff, aes(x = time)) +
  geom_line(aes(y = FVV)) + 
  geom_line(aes(y = prec.per, color = "prec")) +
  geom_line(aes(y = FVV.per * 10 + 150, color="FVV")) +
  ggtitle("Annual Trends, scaled")

ggplot(dff, aes(x = time)) +
 # geom_bar(aes(y = prec / 10, color = "rain"), stat = 'identity') + 
  geom_bar(aes(y = (prec - prec.per) / 10, color = "residuals"), stat = 'identity') +
  geom_line(aes(y = FVV - FVV.per + 1.8)) +
  ggtitle("residual rain and residual VV")